arXiv:1507.01508vl [astro-ph.SR] 6Jul2015 


Astronomy & Astrophysics manuscript no. gaurat_2015 ©ESO 2015 

July 7, 2015 


Evolution of a magnetic field in a differentially rotating radiative 

zone 

M. Gaurat 12 , L. Jouve 1,2 , F. Lignieres 1 - 2 , and T. Gastine 3 

1 Universite de Toulouse, UPS-OMP, Institut de Recherche en Astrophysique et Planetologie, 31028 Toulouse Cedex 4, France 
email: [mathieu.gaurat;laurene.jouve;francois.lignieres]@irap.omp.eu 

2 CNRS, Institut de Recherche en Astrophysique et Planetologie, 14 avenue Edouard Belin, 31400 Toulouse, France 

3 Max-Planck-Institut fiir Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Gottingen, Germany 
e-mail:gas tine® mps.mpg.de 

Received / Accepted 


ABSTRACT 

Context. Recent spectropolarimetric surveys of main-sequence intermediate-mass stars have exhibited a dichotomy in the distribution 
of the observed magnetic field between the kG dipoles of Ap/Bp stars and the sub-Gauss magnetism of Vega and Sirius. 

Aims. We would like to test whether this dichotomy is linked to the stability versus instability of large-scale magnetic configurations 
in differentially rotating radiative zones. 

Methods. We computed the axisymmetric magnetic field obtained from the evolution of a dipolar field threading a differentially 
rotating shell. A full parameter study including various density profiles and initial and boundary conditions was performed with a 
2D numerical code. We then focused on the ratio between the toroidal and poloidal components of the magnetic field and discuss the 
stability of the configurations dominated by the toroidal component using local stability criteria and insights from recent 3D numerical 
simulations. 

Results. The numerical results and a simple model show that the ratio between the toroidal and the poloidal magnetic fields is highest 
after an Alfven crossing time of the initial poloidal field. For high density contrasts, this ratio converges towards an asymptotic 
value that can thus be extrapolated to realistic stellar cases. We then consider the stability of the magnetic configurations to non- 
axisymmetric perturbations and find that configurations dominated by the toroidal component are likely to be unstable if the shear 
strength is significantly higher than the poloidal Alfven frequency. An expression for the critical poloidal field below which magnetic 
fields are likely to be unstable is found and is compared to the lower bound of Ap/Bp magnetic fields. 
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1. Introduction 

In the past decades, the study of the interplay between mag¬ 
netic fields and differential rotation in radiative zones has mainly 
been driven by constraints on stellar rotation. This interaction 
has in particular been invoked to explain the flat rotation pro¬ 
file of the solar radiative zone (see for example Mestel & Weiss 
1987; Charbonneau & MacGregor 1993) or more recently, the 
slow rotation of the core of subgiants and giants revealed by as- 
teroseismology (Deheuvels et al. 2014; Cantiello et al. 2014). 

The evolution of the angular momentum in the presence of 
a magnetic field has mainly been investigated assuming axisym- 
metry and neglecting meridional circulation and turbulence. This 
was the case of the numerical simulations by Charbonneau & 
MacGregor (1992) that considered the spin-up of a radiative stel¬ 
lar zone threaded by an initial poloidal field. The authors pointed 
out the role of the field geometry, more specifically, the exis¬ 
tence of field lines anchored in the core, in reaching or failing 
to reach a solid-body rotation as a final state. They determined 
the timescale necessary to reach this state. This timescale is con¬ 
trolled by the damping of the Alfven waves that are excited by 
the back-reaction of the Lorentz force that follows the winding- 
up of the poloidal field by the differential rotation. Charbonneau 
& MacGregor (1993), Rudiger & Kitchatinov (1996), and Spada 
et al. (2010) studied the spin-down of the radiative zone of the 


Sun under the same physical assumptions. The effect of merid¬ 
ional flows has been addressed in Mestel et al. (1988), who re¬ 
stricted their analysis to the axisymmetric case. The impact of 
non-axisymmetric field components was then studied by Moss 
et al. (1990), Moss (1992) and more recently by Wei & Goodman 
(2015). It was shown that a state of solid-body rotation is reached 
on a timescale shorter than the global diffusive time, regardless 
of the existence of field lines anchored in the core. However, if 
the field strength is small compared to the differential rotation 
strength, the misaligned magnetic field can be axisymmetrised 
before solid-body rotation is reached. 

These studies did not include the effect of possible non- 
axisymmetric instabilities of the magnetic configurations. How¬ 
ever, on the one hand, such instabilities are expected if the 
winding-up of the poloidal field by differential rotation produces 
magnetic configurations that are dominated by the toroidal com¬ 
ponent. Purely toroidal field configurations can indeed be unsta¬ 
ble to various kinds of instabilities, such as the Tayler instability, 
the magneto-rotational instability (MRI), or the buoyancy insta¬ 
bility (see for a review Spruit 1999). On the other hand, magnetic 
fields with toroidal and poloidal components of the same order 
of magnitude can remain stable, as was found in recent numeri¬ 
cal simulations (e.g. Braithwaite & Nordlund 2006; Braithwaite 
2007). For the angular momentum transport, the occurrence of 
such instabilities is a crucial issue since the development of a 
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non-axisymmetric instability could profoundly modify the redis¬ 
tribution of angular momentum. Estimates of the transport re¬ 
sulting from the Tayler instability have been proposed in Spruit 
(2002) on phenomenological grounds, while recent numerical 
simulations by Rudiger et al. (2015) quantified the transport in¬ 
duced by the so-called azimuthal-MRl in the simplified setup of 
a fluid of constant density confined between two rotating cylin¬ 
ders. 

Beyond the problem of angular momentum transport, non- 
axisymmetric instabilities could also affect the value of the mag¬ 
netic field measured at the surface of stars. Notably, spectropo- 
larimetry is an observational technique that allows measuring 
the line-of-sight component of the magnetic field vector inte¬ 
grated over the whole visible surface. This averaged field could 
be strongly decreased as one goes from a well-organized large 
scale magnetic field to a destabilized configuration where oppo¬ 
site polarities cancel each other, especially if the typical length 
scale left by the instability is much smaller than the stellar ra¬ 
dius. This effect has been invoked to explain the magnetic di¬ 
chotomy observed among intermediate-mass stars between the 
strong mostly dipolar Ap/Bp magnetic fields and the ultra-weak 
magnetic fields of other intermediate-mass stars (Auriere et al. 
2007; Lignieres et al. 2014). In this scenario, the observed lower 
bound of Ap/Bp magnetic fields (~ 300 Gauss for the dipolar 
strength) corresponds to the limit between stable and unstable 
magnetic fields. On phenomenological grounds, Auriere et al. 
(2007) approximated the maximum toroidal field produced by 
the winding-up of the poloidal field to — RD. y/4np (where R 
is the radius of the star, £1 is the rotation rate and p is the den¬ 
sity) and assumed that the magnetic configuration is unstable as 
soon as the toroidal field dominates the poloidal field at the stel¬ 
lar surface. They thus derived an order-of-magnitude estimate of 
the critical field B c = RQ yj4np that separates stable and unsta¬ 
ble magnetic fields. This estimate appears to match the observed 
value for a typical Ap/Bp star. 

Our ultimate goal is to investigate the occurrence and the 
effect of non-axisymmetric instabilities of the magnetic field in¬ 
duced by the winding-up of a poloidal field in a differentially 
rotating stellar radiative zone. In the present paper, we use ax- 
isymmetric simulations to explore the different magnetic config¬ 
urations obtained for various differential rotation profiles, den¬ 
sity stratifications, and boundary conditions. In this systematic 
parameter study, we focus on the ratio between the toroidal and 
poloidal magnetic field as a key parameter governing the stabil¬ 
ity of the magnetic configuration. We discuss the occurrence of 
instability in fields that are dominated by their toroidal compo¬ 
nent, for which we use local stability criteria for purely toroidal 
fields as well as results from recent 3D numerical simulations by 
(Jouve et al. 2015) performed for particular configurations. The 
advantage of these 2D axisymmetric numerical simulations over 
3D simulations is to give access to a much broader parameter 
range and in particular to explore asymptotic behaviours in the 
regimes of low magnetic diffusivities and viscosities and large 
density stratifications that are unaccessible to 3D models. 

In Sect. 2, the mathematical formulation of our problem is 
described. We set down the basic assumptions that we adopted 
and then introduce the equations governing the joint evolution of 
the differential rotation and the magnetic field. In Sect. 3, the re¬ 
sults of the numerical simulations are presented. In Sect. 4, sim¬ 
ple models that describe the evolution of the toroidal to poloidal 
field ratio are presented and compared to the numerical results. 
In Sect. 5, we discuss the stability of the magnetic configurations 
obtained by the simulations. 


2. Mathematical formulation 

2.1. Simplifying assumptions and governing equations 

We consider the axisymmetric evolution of an initially poloidal 
field submitted to differential rotation in a spherical shell. Be¬ 
yond the assumption that the flow remains axisymmetric over 
time, we also neglect the meridional circulation. The dynamics is 
therefore not affected by buoyancy effects and is only governed 
by the induction and the azimuthal momentum equations. The 
stellar equilibrium structure comes into play through the den¬ 
sity stratification. We also consider timescales shorter than the 
magnetic diffusion time on a radius length scale. Under the as¬ 
sumptions of axisymmetry and without meridional circulation, 
the initial poloidal field only evolves through Ohmic dissipation 
and is thus considered as constant in time. 

The governing equations are thus 

dB' I 1 \ 

-^ = rsin0(B p .V)Q + i7 A-—-j- B,, , (1) 

at \ r 2 sin-67 

p r sin 0 ^ (B p • V)(r sin6»B^) 

ot Anr sin 0 

+ ^( A ~ 9 1 2 J ^ sin ^) ’ ( 2 ) 

\ r 2 sin- 0) 

where the rotation rate £2(r, 6, t) is related to the azimuthal veloc¬ 
ity by Utpfr, 6,t) — r sin 6 Fl(r, 0, t ) and the magnetic field reads 

B(r, 6, t ) = Bp(r, 6) + B v (r, 0, t) , (3) 

B^fr, 0, t) being its azimuthal component and B p (r, 6), the time- 
independent poloidal component. In these equations, the mag¬ 
netic diffusivity r/ and the dynamic viscosity p are uniform, while 
the density p is given by a polytropic model of a star with a poly¬ 
tropic index n — 3 (Rieutord et al. 2005). 

2.2. From three to two dimensionless numbers 

A possible and straightforward nondimensionalization of these 
equations is obtained by taking the stellar radius R as the refer¬ 
ence length scale, 1/Qo as the reference timescale (where Qq is 
the surface rotation rate at the equator), B ft the surface poloidal 
magnetic field taken at the equator as a magnetic field unit, and 
po the surface density as a density unit. The dynamical Eqs. 
(1) and (2) are then controlled by three different dimensionless 
numbers, namely the Reynolds number Re = £h, R 2 /vo, where 
vo = p/po is the kinematic viscosity at the surface, the Elsasser 
number A = Br./(4npo Q () ?/), and the magnetic Prandtl number 
P m = vo/p. However, it is also possible to reduce the number 
of dimensionless parameters to two by instead employing the 
Alfven timescale t& p = R \[4nfH\/Bo as the reference timescale 
and R Q (l \j4np {) as a magnetic field unit for the toroidal com¬ 
ponent. In addition to the reduction of the number of dimension¬ 
less parameters allowed by this choice, it is motivated by the fact 
that B\p is generated by the Q-effect and hence related to Q () in 
a differentially rotating radiative zone. Bq remains the poloidal 
magnetic field unit and Qo the angular rate unit for this nondi¬ 
mensionalization. The ratio of the toroidal to poloidal magnetic 
fields is then expressed as 

B^ _ B^ tAp^ 

B P 1T P ki 
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where t& — 1 /Qq, and where tildes indicate dimensionless quan¬ 
tities. With this nondimensionalization, Eqs. (1) and (2) read 


<9^, _ - 1 (~ 1 \~ 

-^ = rsin0(B p -V)Q+- A - — \B V , (5) 

at L u \ r 2 sin" 8 / 

_ dh 1 ~ - ~ 

p r sinf? = —-(B p ■ V)(rsin0B„) 

dt r sin (9 


+ 7 ^ (a- 1 2 WsingQ) . ( 6 ) 

L„ \ T^sin 0) 

L u , the Lundquist number, is defined as follows: 

j _ _ R B{) 

L U - - - _ 

l Ap tj 

This reduction from three to two dimensionless numbers greatly 
facilitates the exploration of the parameter space. 



2.3. Initial and boundary conditions 


r^j 


n 


r\j 


n 




Fig. 1: Meridional cut of ff for the radial profile defined by Eq. 13 
(left) and for the cylindrical profile defined by Eq. 15 (right). The dotted 
black lines represent the isocontours of . The white lines represent the 
poloidal magnetic field lines and the red line the one on which Alfven 
waves propagate the slowest. The radius at the inner boundary is r = 
R C = 0.3R 


The initial magnetic field is a dipole (see white lines in Fig. 


1 ): 

2 cos 9 Bo R 3 sin 8 Bo R 3 
B p (r, 6, t = 0) =- r- e r +- t - e g 

7° 7 ° 

BJr, 6,t = 0) = 0 . 


( 8 ) 

(9) 


As stated above, the poloidal component is time-independent, 
while the toroidal component will evolve in time. The boundary 
conditions are a perfect conductor at the bottom boundary and 
an insulating medium at the outer radius. This translates into 


d(rBJr , 9, t)) 


dr 


= 0 and BJr — R,9,f) — 0 


( 10 ) 


r=R c .S.t 


where R, is the radius of the inner boundary. We here explore 
different initial and boundary conditions for the rotation. Three 


different radial differential rotations and two cylindrical profiles 
were considered: 


Q.(r) = Q () 


R 


( 11 ) 


/?2 

Q(r) = fi () — 
r A 


Cl(r) - Qq 


1 - Cl (r - Rc) 2 R/r 3 - C2 (r - R c ) 2 /r 2 

1 -(ci +c 2 )(l -RJR) 2 


( 12 ) 

(13) 


Q(uj) = Qo 

Q(ut) = Qo 


1 + (m/R) 4 


12 


l + ll(m/R) 4 ’ 


(14) 

(15) 


where m = r sin 8 is the distance from the rotation axis. Two 
types of boundary conditions are used at r = R c . namely a stress- 
free or a fixed rotation, while a stress-free boundary condition is 
always used at r — R. These different cases are summarized in 
Tables 3 and 4 and illustrated in Fig. 1 which displays a merid¬ 
ional cut of a radial (left frame) and a cylindrical (right frame) 
differential rotation profile. 

From our reference density profile corresponding to an n — 3 
poly trope, we create different profiles with ratios pjpo ranging 
from 1 to 10 6 (where p c is the density at the radius r — R c ) by 
cutting the polytropic profile at different radii. The created pro¬ 
files are then expanded on the whole computational domain. Fig¬ 
ure 2 shows the radial dependence of the dimensionless Alfven 
velocity va p = B p / yfp for the different density profiles used in 
our simulations. The corresponding time for an Alfven wave to 
propagate from r = R c to r — R along the equator would vary 
from 0.25 t^ P for pjpo - 1 to 22 t Ap for pjpo = 10 6 . 



Fig. 2: vjp at the equator as a function ofTfor different density contrasts 
Pc/po- An = 3 polytropic profile was cut at different radii to produce 
these various pjpo- 
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Fig. 3: The reference case: Meridional cut of B v , B v /B p and Q for different times. t nmx , the time at which B^/Bp is maximal, is equal to 1.85 t^ p 
in this simulation. The colourbar at the bottom of each panel indicates the minimal and maximal values of each variable. The parameters of this 
simulation are L„ = 10 2 , a density contrast of pjpa = 7 x 10 3 , an initial radial differential rotation given by Eq. (13) and a fixed rotation at the 
inner radius. 


2.4. Numerical model 

To solve Eqs. (5) and (6) with the boundary and initial 
conditions, we used the two-dimensional axisymmetric code 
STELEM (Charbonneau & MacGregor 1992; Jouve et al. 2008). 
This code uses a finite-element method in space and a third-order 
Runge-Kutta scheme in time. The computational domain is lim¬ 
ited to the annulus r e [R C ,R], 6 e [0, n], where R, = 0.3 R 
was chosen for all our simulations. P m was fixed to 1 in all our 
simulations while the Lundquist number L u is varied from 10 to 
10 5 . Depending on the value of L u and on the boundary and ini¬ 


tial conditions, different spatial and temporal resolutions are re¬ 
quired. We used a spatial mesh ranging from N r xNg = 128 x 128 
to 512 x 512 nodes depending on the cases. The integration time 
was adjusted to satisfy the Courant-Friedrichs-Lewy criterion. It 
varies from dt = 10 -3 t\ p to dt = 1CF 8 t^ P - 

With the initial conditions used in this study, our problem has 
an equatorial symmetry. This is clearly visible in the meridional 
cut of Fig. 1 for example. However, for comparison to future 
studies with more complex symmetry, we chose to include both 
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hemispheres in our computations and always represent the full 
meridional plane. 

We note that we do not present any results for L u > 245 
with a stress-free boundary condition at r — R c . In fact, with this 
boundary condition there is a Hartmann boundary layer at the 
inner radius. The thickness 5 of this boundary layer is known to 
be (Dormy et al. 1998) 


R Bp 

B p • n L u 


(16) 


where n is the unit vector normal to the boundary. The number 
of mesh points is increased according to Eq. (16) to ensure that 
the Hartmann layer is much larger than the size of the numerical 
grid. However, for large L„, this imposes prohibitive restrictions 
on the time step. 


shift between waves on neighbouring magnetic field lines. As we 
can see in panel e, strong gradients of B p and O are then created 
between these waves. The effect of magnetic and viscous dif¬ 
fusion is thus increased. This is the phase-mixing mechanism as 
described in Ionson (1978), Heyvaerts & Priest (1983), or Parker 
(1991). The oscillations of li f then decrease in amplitude, as 
shown in panel /. After this dissipative phase controlled by the 
phase-mixing mechanism, the system evolves towards a steady 
state known as Ferraro’s law of isorotation (Ferraro 1937), in 
which Q is constant along poloidal field lines. In our case, the 
only possible steady state compatible with the boundary condi¬ 
tions and Ferraro’s law is uniform rotation. This state of uniform 
rotation is reached on timescales of the order of a few diffusive 
timescales which are much longer than the timescales we are in¬ 
terested in here. It is thus not shown in Fig. 3. 


3. Description of the results 

The evolution of the magnetic field and the rotation rate is first 
qualitatively described pointing out the successive phases and 
the associated physical mechanisms. We then focus on the ra¬ 
tio between the toroidal and the poloidal magnetic field as a key 
parameter in the study of the stability of a magnetic field. The 
influence of the Fundquist number, the density contrast, the ini¬ 
tial and the boundary conditions on the highest value reached by 
this ratio is investigated in detail. 

3 .1 . Evolution of the magnetic field and the differential 
rotation in a typical case 

Figure 3 shows the evolution of the toroidal magnetic field B ip , 
the ratio between B p and the poloidal magnetic field B p , and the 
rotation rate Q for a typical case (hereafter called the reference 
case). The parameters of this simulation are L u = 10 2 , a density 
contrast p c /po — 7 X 10 3 , an initial radial differential rotation 
given by Eq. (13) and a fixed rotation at the inner radius. Differ¬ 
ent phases in the evolution of B v and Q can be distinguished. The 
first is the winding-up phase, where the shearing of the poloidal 
field by the differential rotation, the so-called O-effect, generates 
B p without any back-reaction on the differential rotation (panel 
a). Since B p is generated by the term B p • VO, the fact that B. p 
is antisymmetric with respect to the equator directly relates to 
the properties of symmetry of Q and B p . After this phase, the 
Forentz force back-reacts on the differential rotation, leading to 
the propagation of Alfven waves along the poloidal magnetic 
field lines in both directions. In this simulation, we mainly ob¬ 
serve an outward propagation from the internal radius because 
the initial perturbation of the system is concentrated close to the 
bottom boundary. 

In panel d, a region of opposite sign is visible on B. f close to 
the core. It is due to the reflection of Alfven waves on the equa¬ 
tor, where conditions on B p and Q allow for such reflections. 
Indeed, the ability for a boundary to reflect Alfven waves de¬ 
pends on the characteristics of B, f and Q at this boundary. In our 
simulations, the equator behaves as a boundary where B^ — 0 
(insulating condition) and dEl/dd = 0 (stress-free condition). 
This pair of boundary conditions induces a perfect reflection of 
Alfven waves at the equator (Schaeffer et al. 2012). The subse¬ 
quent evolution is characterized by the propagation and the re¬ 
flection of the waves either at the equator, at the surface, or at the 
bottom boundary. 

The variable Alfven velocity and the different distances that 
the Alfven waves have to travel before reflexion lead to a phase 



Fig. 4: Evolution of the spatial maximal value of B p /B p for the ref¬ 
erence case. The horizontal dashed line indicates max{B p lB p ) and the 
vertical dashed line indicates t max . 


3.2. Evolution of the ratio between the toroidal and the 
poloidal magnetic field 

Since non-axisymmetric instabilities are expected to occur when 
the magnetic configuration is dominated by the toroidal com¬ 
ponent, we now focus on the highest value of B^/Bp achieved 
during the winding-up process. In the cases we considered, this 
maximum is always reached before the phase-mixing episode 
takes place (panel c for the simulation shown in Fig. 3). Figure 
4 shows the evolution of the spatial maximal value of B^/Bp for 
the reference case. This quantity is defined as the highest value 
reached in the whole domain. At first, we observe a nearly linear 
growth of the spatial maximal value of B^/Bp. As can be seen 
in panels a, b , and c of Fig. 3, during this phase the location 
of the maximum moves as an Alfven wave along a poloidal field 
line. The B^/Bp ratio reaches its highest value in time and space, 
max(Bp/Bp), at a time denoted t max . At later times in Fig. 4, the 
spatial maximal value of B^/B p oscillates due to reflections of 
Alfven waves at the boundaries and rapidly decreases towards 0 
as the system evolves towards a steady state through the dissipa¬ 
tion induced by the phase-mixing mechanism. The temporal evo¬ 
lution presented in Fig. 4 is typical for most of our simulations. 
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Fig. 5: B v /B p (top panels) and Q (bottom panels) as a function of s/L for different times, s is the curvilinear coordinate of the poloidal field line 
in which max(B f /B p ) is reached, and L the length of this field line, s = 0 locates the inner radius and s = L the equator. The simulation is the same 
as in Fig. 3 with the exception of L„ = 10 4 . The vertical dashed line indicates s max the location of max{B p l B p ) on the poloidal field line. 


We note that (r max , 6 max ), the location at which max(B v /B p ) is 
reached, is visible in panel c of Fig. 3 at mid-radius and low lat¬ 
itude. In this zone, the poloidal magnetic field lines do not reach 
the external radius, they are closed within the shell (see Fig. 1). 

To illustrate the physical mechanisms that give rise to 
max(Bp/B p ), one can look at the evolution of B p / B p and O along 
the poloidal field line on which max{Bp/B p ) is reached. Figure 5 
displays this evolution for the parameters of the reference case, 
except that L u — 10 4 instead of 10 2 . Top panels show the growth 
and decrease of Bp/B p between 0 and 2 t max , max(B lfi /B p ) be¬ 
ing reached at time t max and location s, nnx . The bottom panels 
show that O behaves like a standing wave rather than a trav¬ 
eling wave. This is due to the boundary conditions used in the 
simulation and the spatial scale of the initial differential rota¬ 
tion, which is comparable to the length of the field line. then 
also behaves as a standing wave, but this is not clearly visible in 
the top panels where B ip lB p is represented. Indeed, since B p is 
not uniform along the field line, the standing wave character of 
B^ is not obvious in B p / B p . We also note that at approximately 
t = t max , the gradient of Cl projected onto the poloidal field line, 
dCl/ds = e s ■ VO where ,v is the curvilinear coordinate on the 
field line, changes sign. According to the induction Eq. (5) in 
which the effects of magnetic diffusion are neglected ( L u » 1), 
d(Bp/B p )/dt changes sign together with dCl/ds, and we indeed 
observe in Fig. 5 the decrease of Bp/B p after t max . 

3.3. Summary of the parameter study 

The quantity max{Bp/B p ) defined in the previous section deter¬ 
mines whether the magnetic configuration is locally dominated 
by the toroidal or the poloidal field, and t max gives the time at 
which max(Bp/B p ) is reached. Tables 1, 2, 3, and 4 summa¬ 
rize the values of max(Bp/B p ) and t max obtained by varying the 
Lundquist number L„, the density contrast p c /po, the initial dif¬ 
ferential rotation profile, and the boundary condition for Cl at 
r = R c , respectively. In the following, these results are discussed, 
starting with the effect of L u . The two last columns show the pre¬ 
dictions of models that are presented in the next section. 


3.4. Influence of the diffusivities 




Fig. 6: max(B v /B p ) (top) and t max (bottom) as a function of L„ with the 
physical ingredients used for the reference case. 


Figure 6 displays max(Bp/B p ) and t max as a function of L u 
for the data listed in Table 1. Both max(B f /Bp ) and t max increase 
with L„. We expect maxi B p / B p ) to increase when the effects of 
diffusion on the toroidal magnetic field decrease and thus when 
L u increases. To explain why t max increases with L u , we note 
that the diffusion term in the induction equation is expected to 
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L u 

tmax 

(simulation) 

tmax 

(SW model) 

max(B v / B p ) 
(simulation) 

maxIB^/Bp) 
(Auriere model) 

max(B p /B p ) 
(SW model) 

10 1 

0.76 

2.55 

0.75 

4.86 

4.16 

10 2 

1.85 

2.55 

3.04 

4.86 

4.16 

245 

2.22 

2.55 

4.29 

4.86 

4.16 

10 3 

2.35 

2.55 

4.98 

4.86 

4.16 

10 4 

2.48 

2.55 

5.50 

4.86 

4.16 

10 5 

2.53 

2.55 

5.64 

4.86 

4.16 


Table 1: Influence of L u on t max and max(B (f / B p ) obtained by the simulations and predicted by the model of Auriere 
et al. (2007) and the model of standing waves (SW) presented in Sect. 4.2. The parameters of the simulations are 
Pc! Pa - 7 x 10 3 , the initial radial differential rotation given by Eq. (13) and fi fixed at the inner radius. 


Density contrast 

Pc/Po 

tmax 

(simulation) 

tmax 

(SW model) 

max(B^/Bp) 

(simulation) 

max(B v /B p ) 
(Auriere model) 

maxIB^/Bp) 
(SW model) 

1 

0.57 

0.45 

1.55 

1.00 

1.80 

5.6 x 10 2 

1.15 

1.11 

2.68 

1.63 

1.94 

7 x 10 3 

2.52 

2.55 

5.64 

4.86 

4.16 

2 x 10 4 

3.89 

3.92 

8.55 

7.95 

6.22 

3.5 x 10 4 

4.95 

5.00 

10.8 

10.4 

8.06 

7 x 10 4 

6.67 

6.77 

14.5 

14.4 

11.1 

10 5 

7.86 

7.98 

17.0 

17.2 

12.9 

3 x 10 5 

13.0 

13.3 

28.0 

29.4 

21.7 

10 6 

22.9 

23.5 

49.0 

53.2 

37.7 


Table 2: Influence of the density contrast pci Pa on t„ mx and max{B^/B p ). The parameters of the simulations are L u = 10 5 , the initial 
radial differential rotation given by Eq. (13) and fl fixed at the inner radius. 


Initial condition 

for Q 

tmax 

(simulation) 

tmax 

(SW model) 

max(B^/ Bp) 

(simulation) 

max(B v /B p ) 
(Auriere model) 

max(By/ Bp) 
(SW model) 

Eq. 11 (radial profile) 

2.48 

2.55 

3.63 

4.58 

2.93 

Eq. 12 (radial profile) 

2.50 

2.55 

18.5 

8.85 

13.1 

Eq. 13 (radial profile) 

2.48 

2.55 

5.50 

4.86 

4.16 

Eq. 14 (cylindrical profile) 

0.96 

2.55 

0.31 

3.82 

0.39 

Eq. 15 (cylindrical profile) 

1.35 

2.55 

2.19 

6.13 

2.91 


Table 3: Influence of the initial differential rotation on t max and max(B lf: /B p ). The parameters of the simulations are L„ = 10 4 , 
Pc!Pa - 7 x 10 3 and Q fixed at the inner radius. 


be negative (positive) near local maxima (minima) of B lfi . Thus, 
diffusion contributes to a faster change of sign of dB^/dt and 
hence to a smaller t max . 

More interestingly, we see in this figure that an asymptotic 
value is reached for L u > 10 4 . The same asymptotic behaviour 
with L„ is observed for all cases considered. The fact that the 
magnetic and viscous diffusions do not affect the maximum of 
B^/Bp above a certain L u is not surprising since this maximum 
is reached on a timescale of the order of an Alfven time, this 
timescale becoming much shorter than the diffusive timescale as 
L u increases. This asymptotic behaviour constitutes an interest¬ 
ing feature when applications to realistic stellar radiative zones 
are considered. 

3.5. Influence of the density profile 

The dependence on the density contrast was studied numerically 
by varying pdpo between 1 and 10 6 (see Table 2). For a given 


configuration (L u = 10 5 , an initial differential rotation given by 
Eq. 13 and Q fixed at the inner radius), the two top frames of 
Fig. 7 show that cos 6 max '&n<\T max respectively increases and de¬ 
creases towards an asymptotic value as p c /po increases. Thus, 
the location of the maximum B^/Bp remains unchanged above 
a certain p c /po- As B P is produced through the propagation of 
Aflven waves along field lines, the field line where maxiB^/ B p ) 
is reached is the one that goes through 9 max and r max . This field 
line will remain the same above a certain pdpo ratio. The same 
asymptotic behaviour with pdpo is observed for the other simu¬ 
lations not listed in Table 2, performed for different L„ and initial 
differential rotations. For all the cases considered, the asymptotic 
value of T max is comprised between 0.6 and 0.8, while cos 6 max 
is comprised between 0.3 and 0.4. Hence, max(B ifi / B p ) is always 
reached in a zone located between the middle and the top of the 
radiative zone, close to the equator. 

Table 2 shows that maxiB^jB p ) and t max monotonically in¬ 
crease with pdpo • If we consider max(Bipl B p ) x ^Jpo/p c and 
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Initial condition 

for n 

Boundary condition 
for Cl at r = R, 

tmax 

(simulation) 

tmax 

(SW model) 

max(B p l B p ) 
(simulation) 

max(B f / B p ) 
(Auriere model) 

max(B lf / B p ) 

(SW model) 

Eq. 13 (radial profile) 

fixed value 

2.22 

2.55 

4.29 

4.86 

4.16 

Eq. 13 (radial profile) 

stress-free 

1.69 

1.27 

2.97 

4.86 

3.03 

Eq. 14 (cylindrical profile) 

fixed value 

0.85 

2.55 

0.26 

3.82 

0.39 

Eq. 14 (cylindrical profile) 

stress-free 

0.85 

1.27 

0.26 

3.82 

0.28 


Table 4: Influence of the boundary condition of fi at the inner radius on t max and max(B v /B p ). The parameters of the simulations 
are L„ = 245 and p c /po = 7 x 10 3 . 
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Fig. 7: cos 6 max , r max (two top frames), max(B v /B p )X yp c /p 0 and t max x 
yp c I Pa (two bottom frames) as a function of pdpo- The parameters of 
the simulations are L„ = 10 5 , a prohle of fl defined by Eq. (13) and fi 
fixed at the inner radius. 


tmax X yPol Pc* we find that an asymptotic limit is reached for 
large pdpo. This is shown in the two bottom frames of Fig. 7. In 


other words. 


ma.\(B v /B p ) 

R n 0 yAnpJBa 


Cl and 


R yAnpd Bo 


C2 as pdpo 


tends towards realistic stellar values, where Cl and C2 are inde¬ 
pendent of the density contrast. For the different simulations per¬ 
formed, the value of Cl varies between 3x ICC 3 and 2x 10 1 and 
C2 between 10 2 and 3 x 10 2 . The variations of these constants 
are mainly due to the initial differential rotation considered and, 
to a smaller extent, to the boundary condition adopted at the in¬ 
ner radius. This asymptotic behaviour is expected if the physical 
mechanism that gives rise to max(B if /B p ) is confined in the stel¬ 
lar interior since above a certain p c /po ratio the density profile is 
no longer modified in the internal layers. Indeed, maxiB^/B p ) is 
reached along a field line that remains confined in these internal 
layers. 


3.6. Influence of the initial and boundary conditions 

The initial differential rotation can have a strong effect on 
max(Bp/B p ), but it only slightly modifies t max . Indeed, as 
observed in Table 3, the value of max(B f /B p ) varies from 
18.5 t Ap /tn for the 1/r 2 profile (Eq. 12), which is the strongest 
shear we can use that is hydrodynamically stable, to 0.31 t Ap /to 
for the cylindrical profile given by Eq. (14). Generally, the 
stronger the shear, the larger max(B^,/B p ). But we also find that 


the cylindrical differential rotation induces less B p than the ra¬ 
dial one because the isocontours of O are more aligned with the 
poloidal magnetic field lines in the cylindrical case, causing the 
Q-cffcct to be less efficient. The value of t max is approximately 
independent of the radial differential rotation profile. This indi¬ 
cates that the time to reverse the shear does not seem to depend 
on its intensity. However, t max varies from 0.96 t Ap to 1.35 t Ap for 
the two cylindrical profiles considered. 

Our variables of interest maxiB^/B p ) and t„ mx are less in¬ 
fluenced by the boundary conditions on Q than by all the other 
parameters considered in this study. As illustrated in Table 4, for 
radial differential rotation profiles, mux(B p j B p ) and t max are al¬ 
ways slightly larger for a fixed value than for a stress-free bound¬ 
ary condition, 1.01 to 1.66 and 1.03 to 1.37 times higher, re¬ 
spectively, depending on the cases. For the cylindrical profiles, 
mux(B p )B p ) and t max reach the same values, independently of 
the boundary conditions. 

The boundary conditions other than the one on Q at r = R, 
only have a limited impact on max(B p /B p ) and t max . We thus do 
not show their influence here. If longer timescales were consid¬ 
ered (and final steady states), they would play a significant role, 
however. 


4. Simple models 

In this section two models are presented that provide an estimate 
of maxiB^/B p ) and t max . Their predictions are compared with the 
results of the simulations. The first model is a local model that 
uses the same phenomenological assumptions as Auriere et al. 
(2007). The second model consists of estimating the poloidal 
field line on which max(B p /B p ) is reached and then of solving a 
ID Alfven wave equation along this field line. 


4.1. Local model of Auriere et al. (2007) 


Starting from the ideal induction equation for the toroidal com¬ 
ponent 

8B X 

~^-=r sin0B p -VQ , (17) 

Auriere et al. (2007) assumed that the toroidal field back-reacts 
on the differential rotation after a time T„ mx = (/v Ap , correspond¬ 
ing to the period of an Alfven wave of wavelength t = |VQ|/Q, 
the length scale of the initial differential rotation. Then, by inte¬ 
grating Eq. (17) from t = 0 to t — T rrmx and assuming that the 
differential rotation is constant up to that time, we obtain 


Dtp 

—(r,e,t = T mclx 

Bp 


) = 


r sin0Q(r, 6, t = 0) 
v Ap (r, 6) 


(18) 


Taking the spatial maximum of Btp/B p in Eq. (18), 
maxfBp/Bp), r max , 0 max and then T rnax can be determined. 
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Of the assumptions of this model, neglecting magnetic diffu¬ 
sion appears justified since the numerical results have shown that 
max(By/ B p ) becomes independent of L u when L u > 10 4 - 10 5 . 
Relating T max to the length scale of the initial differential rota¬ 
tion f leads to a simple expression, but then By/B p taken at T rnax 
does not depend on the initial gradient of D and in particular on 
its projection onto poloidal field lines. For example, if Ferraro’s 
law is verified initially, By should remain equal to zero, but the 
model incorrectly predicts that induction takes place. Another 
problem is that with a local model, the effects of the boundary 
conditions are not taken into account. 


4.2. Model of standing Alfven waves on a poloidal magnetic 
field line 

Since the initial perturbations of D and By propagate as Alfven 
waves along B p , a more accurate model would be to solve a one¬ 
dimensional wave equation along the poloidal field line where 
max(Bp/B p ) is reached. As we expect this field line to be the 
one for which the O-effect lasts the longest, we assume that 
maxi By/ B p ) is reached on the field line that maximizes the 
Alfven propagating time over the line length L (taken from the 

inner radius to the equator or the surface) t(L ) = J ds/uA P . 
The field line obtained using this condition for a case with 
Pc/po — 7 X 10 3 is indicated in red in Fig. (1). 

Neglecting the magnetic and viscous diffusions and assum¬ 
ing that the scale of variation of r sin 9 and B p along B p is larger 
than the scale of dD/ds (Mestel & Weiss 1987), the ID wave 
equation to be solved on this field line is 


d 2 D(s, t) 
ds 2 


1 d 2 D(s,t) 




dt 2 


= 0 


(19) 


As already noticed, Q behaves approximately as a standing 
wave for the initial differential rotation considered here (see Fig. 
5). We thus search for normal-mode solutions f(s) exp (itot) of 
the wave equation. In addition, we consider the fundamental 
mode as an approximate solution of the problem because the ini¬ 
tial perturbation is dominated by its large-scale component. The 
equation for f(s) is 


d-fjs) 
ds 2 


VAp(s) 


m = 0 


( 20 ) 


A simple and explicit solution is obtained assuming that va p 
is a constant equal to the mean Alfven velocity UX// along the 
poloidal field line. We then obtain B p j B p (s, t = t max ) by inte¬ 
grating the induction equation up to t max 


By , „ g(s) In s \ 

— CM = t max ) = ACl - = cos -- , 

Bn v A p \2 LI 


( 21 ) 


for a fixed boundary condition at s = 0 and a stress-free bound¬ 
ary condition at s — L, with t max = = and where AD — D(s = 

VAp 

0, t — 0) - Q(s = L, t — 0) and g(s) = r(s) sind(s). 


By AD g(s ) . . _ 

— (S,t = tmax) = ~Z =■ Sin 7T- 
Bp 2 v Ap 


'(a) • 


( 22 ) 


for a stress-free boundary condition at s = 0 and s - L and with 


tmax ~ 2W P - 


A more accurate solution can be obtained using the Wentzel- 
Kramers-Brillouin (WKB) approximation 


By. 

B, 


<s,t — t max ) — AD 


g(s) 


n t(s ) 


-y Iv Ap (s)va p (L) \2 t(L) 


cos 


(23) 


for a fixed boundary condition at s = 0 and a stress-free bound¬ 
ary condition at s — L, with t max = t(L) and where r(s) = f n —. 

JU \)\p 


By _ _ AD 

,. Os t — trnax) — r. 

Bp 2 


g(s) 


sin n 


y/v Ap (s) v Ap (L) \ t(L) 


t(s) 


(24) 


for a stress-free boundary condition at s = 0 and s - L and with 
t — iff) 

l max 2 * 

We now compare this WKB solution and the model of Au¬ 
riere et al. (2007), hereafter called SW model for "standing 
wave" and Auriere model, with the results of the numerical sim¬ 
ulations. 


4.3. Comparison of the models with the numerical 
simulations 

The values of max(By/B p ) and t max given by the models are com¬ 
pared with the results of the simulations in Tables 1, 2, 3 and 4. 

First, we observe that the Auriere model and the SW model 
reproduce the asymptotic behaviour of max(By/B p ) when p c /po 
increases. Indeed, Table 2 indicates that max(By/B p ) is given 
with a mean relative error of 13% by the Auriere model against 
24% for the SW model. In this particular case, the Auriere model 
agrees better with the simulations than the SW model. However, 
in the other simulations that are not listed in this paper it is gen¬ 
erally the opposite trend. The SW model provides very accurate 
predictions of t max with a mean relative error of only 4%, while 
the Auriere model predicts values of t max several orders of mag¬ 
nitude higher than the simulations (we have thus chosen not to 
list them in the tables). 

As expected, the effect of the initial differential rotation is 
much better reproduced by the SW model than by the Auriere 
model, with a mean relative error on max(By/B p ) of 26% and 
280%, respectively (see Table 3). While the Auriere model does 
not take into account the effects of the boundary conditions, the 
SW model estimates max(By/B p ) with a mean relative error of 
only 16% for the different boundary conditions considered (see 
Table 4). Finally, Fig. 8, for the stress-free BC case, and Fig. 5, 
for the fixed BC, indeed show that the evolution of the pertur¬ 
bation is not far from the standing waves prescribed by the SW 
model. 

We conclude that the SW model provides a useful approxi¬ 
mation for the maximum ratio By/B p in a differentially rotating 
star with a dipolar field, provided that the length scale of the ini¬ 
tial differential rotation is not too small compared to the stellar 
radius. 


5. Towards the instability of the magnetic field 

We now turn to the discussion of possible instabilities of the 
magnetic configurations found in our simulations. Figure 9 
shows the configuration obtained at time t max in a typical sim¬ 
ulation. A 3D rendering of the field lines is shown, coloured by 
the total magnetic intensity. We clearly see in this figure that 
we have a complex structure with mixed poloidal and toroidal 
components, even if the toroidal field dominates at mid-latitudes. 
Moreover, differential rotation is still present in the spherical do¬ 
main at this time of the simulation and may strongly influence 
the stability conditions of this complex magnetic field. 

One way to study the stability of such a complex magnetic 
configuration without approximation is to perform 3D numeri¬ 
cal simulations where the full set of MHD equations is solved. 
This approach has been followed in a companion paper (Jouve 
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Fig. 8: fi as a function of s/L for different times, s is the curvilinear coordinate of the poloidal field line on which max(B v /B p ) is reached, and L 
the length of this field line, s = 0 locates the inner radius and s = L the equator. The parameters of the simulation are L u = 245, pdpo = 7 x 10 3 , 
the differential rotation profile given by Eq. (13) and a stress-free boundary condition. 


et al. 2015), although the magnetic configurations studied in 3D 
have been limited to cases of uniform density, cylindrical dif¬ 
ferential rotation, and low L u values. Jouve et al. (2015) found 
that a magnetic instability is triggered and destroys the large- 
scale magnetic field if the ratio of the poloidal Alfven time t Ap 
to the rotation timescale to is sufficiently high. In unstable cases, 
an enhanced transport of angular momentum due to the turbu¬ 
lence induced by the instability is found, confirming the results 
of Rudiger et al. (2015) in cylindrical geometry. This fully 3D 
study also served to test the capacity of approximate local cri¬ 
teria to determine the stability of these complex configurations. 
In particular, it was found that the dispersion relation of Ogilvie 
(2007) could predict the nature of instability and provide rea¬ 
sonable estimates of its growth rates. Here, we used these lo¬ 
cal criteria to predict the stability of the configuration obtained 
through our axisymmetric simulations. As compared to Jouve 
et al. (2015), we therefore study the stability of a wider range of 
configurations, including the effects of the density stratification, 
of a radial differential rotation profile, and of stable stratification. 

5. 1. Local stability analysis under simplifying assumptions 

From numerous previous studies, we know that a magnetic field 
is more likely to be unstable if it is dominated by one of its 
components, namely when it is either purely poloidal (Markey 
& Tayler 1973; Flowers & Ruderman 1977) or purely toroidal 
(Tayler 1973; Pitts & Tayler 1985). In the previous sections, we 
have established in which conditions the initial differential ro¬ 
tation produces a configuration dominated by the toroidal field. 
For the stability analysis, we only consider the situations where 
the toroidal magnetic field is so dominant that the poloidal com¬ 
ponent can be neglected. We return to this assumption in the next 
sub-section. 

A major difficulty in performing a linear stability analysis 
of our magnetic configurations is that the axisymmetric back¬ 
ground field evolves in time. Therefore, an important quantity to 
consider here is the ratio between the growth time of possible 
instabilities and the evolution timescale of the background mag¬ 
netic field. From now on, we assume that this ratio is low and 
thus that the background state can be assumed to be steady. This 
assumption is discussed in the next sub-section. 

The local magnitude of the differential rotation is measured 
by q£l, where q = d In Ll/d In nr is the shear parameter and m is 
the cylindrical radius. For the different profiles considered (see 


B 

77,3* 


10 


I 


Fig. 9: 3D view of the magnetic configuration obtained with our 2D 
simulations at t max . It is computed for L„ = 10 2 , a density contrast 
Pc/po = 7 x 10 3 , the radial differential rotation profile defined by Eq. 
(13) and Q fixed at the inner radius. 


Eqs. 11 to 15), the shear parameter is always of order unity, 
meaning that the differential rotation is of the order of the lo¬ 
cal rotation rate. 

Under these assumptions, we may proceed to a local lin¬ 
ear stability analysis, using the dispersion relation derived by 
Ogilvie (2007) for the case of a purely toroidal field under the 
influence of a differential rotation, both possessing arbitrary pro¬ 
files. Two instabilities may be expected in this situation: the 
Tayler instability (TI), which derives from free magnetic en¬ 
ergy, and the magneto-rotational instability (MRI), whose source 
of energy is the free kinetic energy of the differential rotation. 
The steady and axisymmetric basic state consists of a differen¬ 
tial rotation profile Ll(w, z) and a purely toroidal magnetic field 
B^ziT, z), where z is the coordinate along the rotation axis. If we 
consider a perturbation such that the displacement is in the direc- 
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tion e, we obtain the following form of the dispersion relation: 


2 2 
m v A 

or -- (e • e r ) 2 N 2 ■ 

w l 


2z&Qe • VjQ + 


2 2 

2 m % 




2 £2 + 


ZCT 2 


2 c'a if 

yj^np 


e-V 



(25) 


where w is the frequency of the perturbation, in is the azimuthal 
wavenumber, va< p = B^/ yjAnp is the toroidal Alfven velocity, 
and N is the Brunt-Vaisala frequency taken as uniform in the 
whole domain. As noted previously, the buoyancy force does 
not affect the evolution of the axisymmetric magnetic field since 
the meridional flows are neglected. It influences the stability of 
the magnetic configurations through the value of the N/Cl ra¬ 
tio, however. As noted by Jouve et al. (2015), another important 
quantity of the magnetic stability is the ratio of the rotation rate 
£1 to the toroidal Alfven frequency <i) Alf = u Ap l r. 



Fig. 10: Contours of the toroidal magnetic field (white lines) and of the 
growth rate of the instability (coloured contours) for the most unstable 
perturbation (m = 7 for the top frame and m = 1 for the bottom frame) 
which is in the direction 6 = 6^. The parameters of the simulations are 
N = 0 and, respectively for the top and the bottom frame Q.la> Aip =10 
and Cl/a)A V * 1 at the maximum of the toroidal field. 


For a given configuration B^(m, z, t ma x) and Q.(m,z,t ma x), 
corresponding to a £i(w, z„ t max )l u) Aip {m, z, t max ) ratio, obtained 
with the 2D numerical simulations and a value of N/Sl, the dis¬ 
persion relation provides the growth rate cr = 3 (w(tjt, z)) of a 
local perturbation characterized by m and e. Two cases are illus¬ 
trated in Fig. 10. In this figure, only the northern hemisphere is 
shown since the results are symmetric with respect to the equa¬ 
tor. The contours of the toroidal magnetic field are drawn (white 
lines) together with the contours of the growth rate of the in¬ 
stability (coloured contours) for the most unstable perturbation, 
which is in the direction e = The only difference in the initial 
setup of these two simulations is the differential rotation profile 
(Eq.14 was used for the case in the top panel and Eq.13 for the 
bottom panel). As a consequence of a different winding-up be¬ 
tween these two cases, the ratio reaches very different 

values. As shown in Jouve et al. (2015), the nature of the insta¬ 
bility that is likely to be triggered in these magnetic configura¬ 
tions depends on this CI/uia^ frequency ratio. For the case in the 
top panel « 10 at the maximum of toroidal field, the dif¬ 

ferential rotation thus plays a dominant role, and the most vigor¬ 
ous instability is the MRI, the most unstable mode having an az¬ 
imuthal wave number in = 7. In contrast, for the other case where 
~ 1 at the maximum of toroidal field, the current-driven 
Tayler instability is triggered and favours the m = 1 mode. The 
location of the instability is also quite different in both cases. In 
the MRI case, the most unstable region is concentrated at rather 
low latitudes and very extended in radius. In the TI case, how¬ 
ever, the unstable zones are mainly located on strong gradients 
of the field close to the bottom of the domain at latitudes 30° 
and 60°. In all cases studied, the value of D./a) A ip is found to vary 
between about 1 and 10 depending on the initial and boundary 
conditions. We thus expect our axisymmetric configurations to 
be subject to an instability, the nature of which only depends 
on this O /cijAtp ratio. Both the MRI and the TI are found to be 
likely to exist in our situations. The location of the instability, 
the growth rates as well as the typical lengthscale of the most 
unstable mode will then differ and the consequences on the ob¬ 
served surface field may be very different between the TI and the 
MRI cases. 

We showed in Sect. 3.5 that when the density contrast is in¬ 
creased between the top and bottom of our domain, both the 
highest value of B f /B p and its location tend to an asymptotic 
value. For the stability conditions, we find that for the case where 
the in = 1 Tayler instability is favoured (bottom panel of Fig. 10), 
increasing the density contrast does not significantly vary the na¬ 
ture, location, and growth rate of the instability. Indeed, as p c /po 
is increased from 7 x 10 3 (case shown in the figure) to 2 x 10 4 
and to 10 6 , the whole magnetic configuration is left mostly un¬ 
changed and the value of D./oj Alf increases from 1 to 1.25 and to 
1.54. In these three cases, the in = 1 mode is found to be the most 
unstable to the TI and is always located at the base of the domain 
in two distinct regions in latitude (i.e. at 30" and 60° where the 
gradients are the strongest). The highest growth rate is always 
around cr » 7 -9 where £l A <p = B Vm /R ^Anp c and B iPm is the 
maximal toroidal field. This is an interesting feature of this work 
since it implies that our results can be extrapolated to realistic 
values of the density contrast not only for the location and value 
of max(B ip /B p ), but also for the instabilities that are likely to be 
triggered in the stellar interior. 

The effect of the stratification has also been explored by 
varying the parameter N/Q. o, where Q (l is the rotation rate at 
the surface. For the case where the Tayler instability is favoured, 
increasing N/ Q<> from 0 (case shown in the bottom panel of Fig. 
10 ) to 8 and to 16 decreases the maximum growth rate of the in- 
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stability from 9.2 Q A ,r to 8.2 f} Alf and to 5.2 £l Alf . The location of 
the instability remains unchanged. For the case where the MRI 
is favoured, the highest growth rate strongly decreases, the most 
favoured wave number m decreases when N /Do increases, and 
the instability totally disappears for N/Sl o > 2. This agrees with 
Spruit (1999), who stated that the MRI is more affected by the 
stable stratification than the TI. It also agrees with Masada et al. 
(2006), who showed that the most unstable azimuthal wave num¬ 
ber m decreases when N /0 () increases. Nevertheless, with a re¬ 
alistic stable stratification (N/CIq - 10), the most unstable per¬ 
turbations are no longer in the direction e = e m , but rather in the 
direction e = eg. Indeed, such horizontal perturbations are not 
affected by the stratification, and magnetic instabilities are still 
present in these cases. For the same case as the one shown in the 
bottom panel of Fig. 10, the m — 1 TI still possesses a signifi¬ 
cant growth rate cr ^ 9 in the direction e = eg. The MRI can 
still be triggered for e = eg, but with a growth rate cr ^ 0.4 
much lower than in the unstratified cases, where cr =* 2.6 Q Aip 
(case shown in the top panel of Fig. 10). We thus still expect 
non-axisymmetric instabilities to develop in the stably stratified 
cases, the main difference being that when CI/loa^, > 1 (MRI 
regime), the growth rate is significantly reduced. 

5.2. Discussion of the stability condition 

One of the simplifying assumptions in the previous section was 
to consider that the instability would grow on a steady axisym- 
metric background state. However, both the background mag¬ 
netic field and the differential rotation always oscillate on the 
short timescales considered in our study. 

From Fig. 4, this oscillation frequency can be estimated to 
be 1 /2 t max . To freely develop, an instability with a growth rate 
cr should therefore fulfil the condition 

cr > ——— . (26) 

2 tmax 

As discussed in the previous sub-section, if Cl/u> Aip > 1, the 
MRI is favoured. Its maximum growth rate is related to the shear 
parameter and the rotation rate in the following way (see for ex¬ 
ample Balbus & Hawley 1991): 


WMRi-qQ /2 ■ (27) 

If Q./oj Alf < 1, the TI is favoured and the growth rate reads 
ctti = ■ (28) 

From the standing Alfven wave model presented in Sect. 4.2, 
the toroidal Alfven frequency (o Aif and t max are given by Eqs. 
(23) and (24), or by Eqs. (21) and (22) in a simplified case. The 
quantity AQ, which appears in these equations, can be expressed 
in terms of the rotation rate O and the shear parameter q\ AQ ~ 
m dQ. I dm = q Q. The condition for the development of a non- 
axisymmetric instability, Eq. (26), can then be rewritten as 


VAp 

AO > — 


(29) 


which is valid for both the MRI and the TI cases. 

In practice, effects ignored in this estimation will make this 
condition more restrictive. First, as shown in the previous sub¬ 
section, the stable stratification is expected to significantly de¬ 
crease the instability growth rate in the MRI case. Second, to 
derive (29), we assumed that the toroidal component dominates 
the poloidal component. This assumption does not impose an ad¬ 
ditional constraint since, according to the SW model, condition 


(29) also ensures that the toroidal component dominates. Nev¬ 
ertheless, even in the presence of a small poloidal component, a 
more conservative instability condition might be necessary (see 
for the TI case Linton et al. 1996; Braithwaite 2009). Finally, 
the 3D simulations by Jouve et al. (2015) indicate that a condi¬ 
tion cr > 1 0y4— is more appropriate than cr > to account 
for the effect of the oscillating background state. Indeed, in the 
regime < cr < 10 ^ 7 — perturbations do grow exponentially, 
but are killed by the background field reversal before they reach 
the level of energy of the axisymmetric field. Taking these three 
effects into account, the condition for the non-axisymmetric in¬ 
stabilities therefore is rather AQ » ^ In the particular case 
of the 3D simulations by Jouve et al. (2015), the condition is 

AQ > 100^. 

5.3. Discussion of the critical surface field 

When written in terms of quantities observables at the stellar 
surface, the condition for instability reads 

A . (30) 

ROo sj^rrpQ ^0 Va p 

where, as we just discussed, the factor C takes into account ef¬ 
fects on the instability ignored in the simple condition (29). In 
the case studied by 3D numerical simulations (where p is uni¬ 
form, AO/Do ~ 1 and the initial field is dipolar), it is found to 
be ~ 10 2 . This expression defines a surface critical field, B cnt , 
below which non-axisymmetric instabilities can be triggered by 
differential rotation. This field is proportional to ROo \/4vrpo and 
depends on internal properties, such as differential rotation and 
field geometry, through the ratio AQ/Do and u Ap (R)/vAp. The ac¬ 
tual value of AQ/Q ( ) will depend on the mechanism that enforces 
differential rotation in the star and might be much lower than the 
AQ/Q() ~ 10 1 - 1 values we considered here as initial condi¬ 
tions. 

Following Auriere et al. (2007), the magnetic dichotomy of 
intermediate-mass stars would be due to the development of non- 
axisymmetric instabilities separating stable strong field configu¬ 
rations, the Ap/Bp stars, from unstable weaker field configura¬ 
tions whose surface average field becomes very weak after the 
destabilization. Accordingly, the lower limit of Ap/Bp magnetic 
fields would be given by a stability condition such as (30) (the 
present stability condition is global and thus differs from the lo¬ 
cal condition found by Auriere et al. (2007)). As we know that 
the observed lower bound of Ap/Bp dipolar fields is ^ 300 G 
and that RDo \J4-71 fig - 300 G for a typical Ap star (R = 3 R 0 , 
P = 5 days, log g = 4dex and 7/n = 10 4 K), the condition for 
instability (30) would match observations if > 1. A 

density profile calculated for such a star thanks to the stellar evo¬ 
lutionary code MESA yields ^ 2 x 10 4 , which implies that 

AQ/Q 0 > 5 X 10 \ However, in making this estimate, we as¬ 
sumed that any non-axisymmetric instability taking place deep 
inside the star strongly affects the observed surface field. This 
is not necessarily true, in which case higher AQ/Q 0 would be 
compatible with the Auriere scenario. 

Given the present uncertainties on the stability condition (see 
previous sub-section), the level of the differential rotation and 
the effect of an internal instability on the surface field, the quan¬ 
titative comparison between the theoretical critical field and the 
observed lower bound of Ap/Bp fields is still rather approxima- 
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tive. More realistic 3D simulations will help, although the den¬ 
sity contrasts encountered in a stellar envelope remain beyond 
reach for such simulations. It will be more relevant to compare 
the dependence of the critical field on the surface rotation rate 
with observations. Existing data indeed point towards a linear 
dependency of the lower bound of Ap magnetic fields (Lignieres 
et al. 2014) compatible with the stability condition (30). 

6. Conclusion 

We performed an extensive study of the evolution of a magnetic 
field in a differentially rotating radiative zone. Among other re¬ 
sults, we provided a systematic estimate of the ratio between the 
toroidal and the poloidal magnetic field. We tested the influence 
of the diffusivities and of the density profile and characterized 
the impact of different boundary conditions and differential rota¬ 
tion profiles. The large parametric study that 2D numerical sim¬ 
ulations made possible has allowed finding asymptotic regimes 
for the diffusion and the density contrast. This in turn enabled 
us to extrapolate our numerical results to realistic stellar diffu¬ 
sivities and density contrasts. A simple one-dimensional stand¬ 
ing wave model was also found to provide reliable estimates of 
max(B, f /B p ) and the time at which this maximum is reached. 

We then discussed the stability of the magnetic configura¬ 
tions dominated by the toroidal component using a local dis¬ 
persion relation derived by Ogilvie (2007) and results of full 
3D numerical simulations performed in a simplified case (Jouve 
et al. 2015). We argued that these magnetic configurations are 
likely to be subject to a magneto-rotational instability or to 
the Tayler instability, the nature of the instability depending 
solely on the ratio of the rotation rate to the toroidal Alfven 
frequency Q./o>a ■ Such non-axisymmetric instabilities are ex¬ 
pected to modify the large-scale axisymmetric configuration 
when the surface poloidal field is weaker than a threshold value 
B cr (j. In the context of the Auriere scenario to explain the di¬ 
chotomy of intermediate-mass star magnetism, our estimate of 
B cr it is compatible with the lower bound of Ap/Bp magnetic 
fields, although a precise quantitative comparison remains dif¬ 
ficult at this stage. 

We neglected meridional flows in our simulations. The effect 
of a prescribed circulation has been considered by Mestel et al. 
(1988) and Moss et al. (1990). They found that it is negligible 
as long as the circulation velocity is very low compared to the 
Alfven speed, which is the case for typical Eddington-Sweet- 
type circulations. Another hypothesis in this work is that of ax- 
isymmetry, when it is known from observations that Ap/Bp stars 
are generally oblique rotators. Non-axisymmetric components of 
the initial poloidal field could have an influence on the stabil¬ 
ity of the magnetic field. Moss et al. (1990), Moss (1992), and 
Wei & Goodman (2015) have shown that if Bq/RCIq yjAjtpQ «c 1 
the magnetic field is symmetrized before the differential rota¬ 
tion decays. In such extreme cases, the magnetic configurations 
are predicted to be unstable by (30). In the opposite regime, 
Bq/RQ) \j4npo » 1, the strong magnetic field quickly sup¬ 
presses the differential rotation and the magnetic field remains 
inclined with respect to the rotation axis. These situations are 
similar to axisymmetric cases in which magnetic configurations 
are also predicted to be stable. In intermediate cases, we do not 
know the effect of the non-axisymmetric components of the field 
on the stability conditions. This question needs to be addressed. 

To proceed, we need 3D numerical simulations that include 
the effects of the stable stratification as well as the mechanism 
that forces the differential rotation. The pre-main-sequence con¬ 


traction phase can in principle generate the required differential 
rotation (Lignieres et al. 2014). While realistic diffusivities and 
density contrasts are beyond what can be reached by 3D simu¬ 
lations, we expect that the dependency of B cnl with the rotation 
rate will not strongly depend on these parameters. The relation 
between B cnt and the surface rotation will then be compared to 
observational constraints on the lower bound of Ap/Bp magnetic 
fields. Another application of such simulations is to consider the 
role of non-axisymmetric instabilities on the efficiency of the an¬ 
gular momentum transport in stellar radiative zones. 
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